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Abstract 

In a previous article |l] we have shown how one can employ Artificial Neu- 
ral Networks (ANNs) in order to solve non-homogeneous ordinary and partial 
differential equations. In the present work we consider the solution of eigen- 
value problems for differential and integrodifferential operators, using ANNs. 
We start by considering the Schrodinger equation for the Morse potential that 
has an analytically known solution, to test the accuracy of the method. We 
then proceed with the Schrodinger and the Dirac equations for a muonic atom, 
as well as with a non-local Schrodinger integrodifferential equation that mod- 
els the n + a system in the framework of the resonating group method. In two 
dimensions we consider the well studied [0 Henon-Heiles Hamiltonian and in 
three dimensions the model problem of three coupled anharmonic oscillators. 
The method in all of the treated cases proved to be highly accurate, robust 
and efficient. Hence it is a promising tool for tackling problems of higher com- 
plexity and dimensionality. 

PACS'96 codes: 02.60.Lj, 02.60. Nm, 02.70.Jn, 03.65.Ge 
Keywords: neural networks, eigenvalue problems, Schrodinger, Dirac, collo- 
cation, optimization. 
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1 Introduction 



In a previous work [|TJ a general method has been presented for solving both or- 
dinary differential equations (ODEs) and partial differential equations (PDEs). 
This method relies on the function approximation capabilities of feedforward 
neural networks and leads to the construction of a solution written in a dif- 
ferentiable, closed analytic form. The trial solution is suitably written so as 
to satisfy the appropriate initial/boundary conditions and employs a feedfor- 
ward neural network as the main approximation element. The parameters of 
the network (weights and biases) are then adjusted so as to minimize a suit- 
able error function, which in turn is equivalent to satisfying the differential 
equation at selected points in the definition domain. 

There are many results both theoretical and experimental that testify for 
the approximation capabilities of neural networks M, ^ M. The most im- 
portant one is that a feedforward neural network with one hidden layer can 
approximate any function to arbitrary accuracy by appropriately increasing 
the number of units in the hidden layer ||. This fact has led us to consider 
this type of network architecture as a candidate model for treating differential 
equations. In fact the employment of neural networks as a tool for solving 
differential equations has many attractive features JT|: 

• The solution via neural networks is a differentiable, closed analytic form 
easily used in any subsequent calculation with superior interpolation ca- 
pabilities. 

• Compact solution models are obtained due to the small number of re- 
quired parameters. This fact also results in low memory demands. 

• There is the possibility of direct hardware implementation of the method 
on specialized VLSI chips called neuroprocessors. In such there 
will be a tremendous increase in the processing speed that will offer the 
opportunity to tackle many difficult high-dimensional problems requiring 
a large number of grid points. Alternatively, it is also possible for the 
proposed method to be efficiently implemented on parallel architectures. 

In this paper we present a novel technique for solving eigenvalue problems of 
differential and integrodifferential operators, in one, two and three dimensions, 



2 



that is based on the use of MLPs for the parametrization of the solution, 
on the collocation method for the formulation of the error function and on 
optimization procedures. 

All the problems we tackle come from the field of Quantum Mechanics, i.e. we 
solve mainly Schrodinger problems and we have applied the same technique 
to the Dirac equation that is reduced to a system of coupled ODEs. In addi- 
tion, for the Schrodinger equation one can employ the Raleigh-Ritz variational 
principle, where again the variational trial wavefunction is parametrized us- 
ing MLPs. For the two-dimensional Hennon-Heiles potential, we compare the 
resulting variational and the collocation solutions. 

A description of the general formulation of the proposed approach is presented 
in section 2. Section 3 illustrates several cases of problems where the proposed 
technique has been applied along with details concerning the implementation 
of the method and the accuracy of the obtained solution. In addition, in a 
two dimensional problem, we provide a comparison of our results with those 
obtained by a solution based on finite elements. Finally, section 4 contains 
conclusions and directions for future research. 

2 The Method 

Consider the following differential equation: 



where if is a linear differential operator, f(r) is a known function, D C R 3 and 
dD is the boundary of D. Moreover, we denote D = DU dD. We assume that 
/ e C{D) and the solution *(r) belongs to C k (D), the space of continuous 
functions with continuous partial derivatives up to k order inclusive (k is the 
higher order derivative appearing in the operator H, H^(f) e C(D)). The set 
of the admissible functions 



HV(r) = /(f), in D 



(1) 



tf(f) = 0, on dD 



(2) 



{ty(r)eC k (D), feDcR 3 , *(f) = on dD} 
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forms a linear space. In the present analysis we also assume that the domain 
under consideration D is bounded and its boundary dD is sufficiently smooth 
(Lipschitzian). 

In order to solve this problem we have proposed a technique Q, that consid- 
ers a trial solution of the form \&t(f) = A(f) + B(f, X)N(f,p) which employs 
a feedforward neural network with parameter vector p (to be adjusted). The 
parameter vector A should also be adjusted during minimization. The specifi- 
cation of functions A and B should be done so that ty t satisfies the boundary 
conditions regardless of the values of p and A. 

To obtain a solution to the above differential equation the collocation method 
has been employed || which assumes a discretization of the domain D into a 
set points rj. The problem is then transformed into a minimization one with 
respect to the parameter vectors p and A: 



mm 

v 



. x U H Mn) - fin)) 2 (3) 



If the obtained minimum has a value close to zero, then we consider that an 
approximate solution has been recovered. 

Consider now the case of the following general eigenvalue problem: 

H^(r) = £*(f), in D (4) 

#(r) = 0, on dD (5) 

In this CcLSG cL trial solution may take the form: \l/t(r) = B(f, X)N(f,p) where 
B(f, A) is zero on dD, for a range of values of A. By discretizing the domain, 
the problem is transformed to minimizing the following error quantity, with 
respect to the parameters p and A: 

Error{p, A) = T\%W ^ 



where e is computed as: 

J m*Hm t dr 
e ~ J\^ t \ 2 df 



(7) 



A method similar in spirit has been proposed long ago by Frost et al 
and is known as the "Local Energy Method". In the proposed approach the 
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trial solution employs a feedforward neural network and more specifically 
a multilayer perceptron (MLP). The parameter vector p corresponds to the 
weights and biases of the neural architecture. Although it is possible for the 
MLP to have many hidden layers we have considered here the simple case 
of single hidden layer MLPs, which have been proved adequate for our test 
problems. 

Consider a multilayer perceptron with n input units, one hidden layer with 
m sigmoid units and a linear output unit (Figure 1). The extension to the 
case of more than one hidden layers can be obtained accordingly. For a given 
input vector r = (r\, . . . ,r n ) the output of the network is N = YhLi Via(zi) 
where Zi = J2]=i w^Tj+Ui, Wij denotes the weight from the input unit j to the 
hidden unit i, Vi the weight from the hidden unit i to the output, Ui the bias of 
hidden unit % and cr(z) the sigmoid transfer function: a(z) = 1/(1 + exp(— z)). 
It is straightforward to show that [|l|] : 

d k N ™ k (fc ) 

aT-LWi ( 8 ) 
Ul j i=i 

where <7j = cr(zj) and a^ k ' denotes the k th order derivative of the sigmoid. 
Moreover it is readily verifiable that: 



dr Xl dr\ 2 " ' drj 



i=l 

where 

n 

A 



-N = Y,vM K) (9) 



Pi = U <k (io) 



/c=i 



and A = £r=i \. 



Once the derivative of the error with respect to the network parameters has 
been defined it is then straightforward to employ almost any minimization 
technique. For example it is possible to use either the steepest descent (i.e. the 
backpropagation algorithm or any of its variants), or the conjugate gradient 
method or other techniques proposed in the literature. We used the MERLIN 
optimization package B, |9| for our experiments, where many algorithms are 
available. We mention in passing that the BFGS method has demonstrated 
outstanding performance. Note that for a given grid point the calculation of 
the gradient of each network with respect to the adjustable parameters, lends 
itself to parallel computation. 
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Using the above approach it is possible to calculate any number of states. This 
is done by projecting out from the trial wavefunction the already computed 
levels. 

If l^o >j |^i >>••■) |*fc > are computed orthonormal states, a trial state > 
orthogonal to all of them can be obtained by projecting out their components 
from a general function \ty t > that respects the boundary conditions, namely: 

|* t >= (1 - |^o >< *o|)(l - l*i >< *i|) • • • (1 - |*fc >< * fc |)|*t > 

= (1 - |^o >< *o| - 1*1 >< *l| . . . - >< > 



3 Examples 



3.1 Schrddinger equation for the Morse Potential 



The Morse Hamiltonian for the I2 - molecule in the atomic units system, is 
given by: 

where V(x) = D[ e - 2ax -2e- ax + 1] and D = 0.0224, a = 0.9374, fi = 119406. 



The energy levels are known analytically [P^j , and are given by: 

e n = ( n + |)(1 - with C = 156.047612535, £ = 5.741837286 10~ 4 . The 

ground state energy is e = 0.286171979 10~ 3 . We parametrize as: 

(f) t (x) = e~^N(x, u,w,v), (3 > 

with N being a feedforward artificial neural network with one hidden layer 
and m sigmoid hidden units, ie: 

m 

N(x, U, W, V) — ^ VjO~{WjX + Uj) 
3=1 

We solve the problem in the interval —1 < r < 2 using 150 equidistant grid 
points with m — 8. We minimize the quantity: 



J(pf{x)dx 
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where e = - ^j^^fa'** ■ We find for the ground state energy the value 
0.286171981 1CT 3 which is in excelent agreement with the exact analytical 
result. 



3.2 Schrddinger equation for muonic atoms 

The s-state equation for the reduced radial wavefunction <f)(r) = rR(r) of a 
muon in the field of a nucleus is: 

h 2 d 2 ^ 



(r) + V(r)(f)(r) = e(j)(r) 



2p dr 2 

with: 0(r = 0) = and 0(r) ~ e~ kr , k > for a bound state. 

u is the reduced muon mass given by: - = — + ^ — hr- r — , where m u is the 

" ° J fj, m M Zm p +Nm n ' M 

muon mass and m p , m n the masses of the proton and neutron respectively. 
Z is the number of protons and N the number of neutrons for the nucleus 
under consideration. (In our example we calculate the muonic wavefunction 
in 82 P6 208 ). 

The potential has two parts, i.e.: V(r) = V e (r) + V p (r), where 

Ve (r) = -e 2 f -^X-A' 
J j|^ — 

is the electrostatic potential, p(r) is the proton number-density given by 

p(r) = A/(l + e {r ~ b)/c ) 
with A = 0.0614932, b = 6.685 and c = 0.545 and 



V L (r) - \v e {r) 
6 



is the effective potential due to vacuum polarization |TU] with a = 137 037 the 
fine-structure constant. 

V L (r) = -2tt— / p{r')r' {\r - r'\[ln{C\r - r'\\ e ) - 1] - (r + r')[ln{C{r + r')/A e - 1]} dr' 



r Jo 

with C = 1.781 and A e the electron Compton wavelength divided by 2n. 

We parametrized the trial wavefunction as: 

<fi t (r) = re _/3r A^(r, u,w,v), (3 > 
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where again N is again a feedforward artificial neural network with one hidden 
layer having 8 sigmoid hidden units. 

The energy eigenvalue is calculated as: 

The intergrals have been calculated using the Gauss- Legendre rule. We used 
80 points in the range [0,40]. The quantity: 

1 ( h 2 d 2 ) 2 

is being minimized with respect to u, w, v. 

We used for the same points as in the Gauss- Legendre Integration. We 
obtained for the energy e = — 10.47MeV\ The radial wavefunction -<j)(r) is 
shown in Fig. 2c. 
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3.3 Dime equation for muonic atoms 



The relativistic Dirac s-state equations for the small and large parts of the 



reduced radial wavefunction of a muon bound by a nucleus are |1 I 

*f(r) + I/( r ) = l(^c 2 - E + V(r))g(r) 
dr r he 

jL g ( r ) _ i^( r ) = _L (/iC 2 + £ _ I/( r ))/(r) 
with [l and V(r) being as in the previous example. 

The total energy E is calculated by: 
1 



E 



JSVM - /»<r)]*^f ^< r > + ™l* + jT V(r)b»(r)-Ar)]*} 



We parametrized the trial solutions /t(r) and <?t(r) as: 

/t(r) = re _/3r iV(r, u), w f , v)), (3 > 
flf t (r) = re _/3r iV(r, w 9 , w g , v g ), (3 > 
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and minimized the following error quantity: 

J2 r[ d /( r ') _|_ fin) _ MC 2 --E+V(r») ^^ r jj2 _|_ \ dgjn) _ g( ri ) _ i ic 2 +E-V(r i ) y^j2 j 

n 2 (r) + f 2 (r]]dr 

The binding energy is given by e = E — fic 2 . We find e = —10.536 MeV. 
The small and the large parts of the radial wavefunction ~/(r) and -g(r) are 
shown in Fig. 2a and 2b, along with the Schrodinger radial wavefunction (Fig. 
2c). The integrals and the training were performed using the same points as 
in the previous example. 

3.4 Non-Local Schrodinger equation for the n + a system 



We consider here the non-local Schrodinger equation : 

h 2 d 2 d> f°° 

V) + V(r)4>(r) + / K (r,r')4>(r')dr' = e<f>(r) 



2/j, dr 2 Jo 

with V{r) = -Voe~P r \ where V = 41.28386, (3 = 0.2751965 and 

K (r,r') = -Ae~^ r2+r ' 2 \e 2krr ' - e ~ 2krr ') with A = -62.03772 , 7 = -0.8025, 

k = 0.46. This describes the n + a system and is derived in the framework of 



the Resonating Group Method |12l , /1 is the system's reduced mass given by: 



1 = J- + 1 

fi m n 2m„+2m p 

We parametrized the trial wavefunction as: 

4>t{r) = re _/3r A^(r, u, w, v), (3 > 

where the neural architecture is the same as in the previous cases and mini- 
mized the following error quantity: 

gj + V{n)Un) + SS° K (r i y)Mr')dr f - eUn)} 2 

JS°tf(r)dr 
where the energy is estimated by: 

jg/2/i ^{ d -t?dr + £ V(r)<p 2 (r)dr + / °° j °° K (r, r'^^drdr' 

We have considered 100 equidistant points in [0, 12] and the computed ground 
state is depicted in Fig. 3, while the corresponding eigenvalue was found equal 
to -24.07644, in agreement with previous calculations M. 
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3.5 Two dimensional Schrddinger equation 



We consider here the well studied || example of the Henon-Heiles potential. 

The Hamiltonian is written as: 

1 d 2 d 2 
H = -2^ + 8f ) + V{X > y) 
with V{x, y) = \{x 2 + y 2 ) + j^(xy 2 - |a; 3 ). 

We parametrize the trial solution as: 

^(z, y) = e- x{x2+y2) N{x, y, u, w {x \w {y \v), A > 

where N is a feedforward neural network with one hidden layer (with m = 8 
sigmoid hidden units) and two input nodes (accepting the x and y values). 



N(x, y, u, vj( x \ w^ y \ v) = ^2 VjCf{xw ( f S> + yuij + Uj) 



We have considered a grid of 20 x 20 points in [—6, 6] x [—6, 6]. The quantity 
minimized is: 

/oo /*oo 
/ dxdy(f) 2 t (x,y) (11) 
-oo J-co 

where the energy is calculated by: 

= IZo S-oo M x , y) H <t>t{x, y)dxdy 

p 00 JZ 4> 2 t (x,y)dxdy 

For this problem we calculate not only the ground state but a few more levels. 
The way we followed is the extraction from the trial wavefunction of the 
already computed levels as described in Section 2. If for example by (j)o(x,y) 
we denote the normalized ground state, the trial wavefunction to be used for 
the computation of another level would be: 

_ /"OO /"OO _ 

<t>t{x, y) = <p t (x, y) - <f> (x, y) / <j) (x', y')M x 'i y')dx'dy' 

J — OO J —CO 

where (j) t (x,y) is parametrized in the same way as before. 

Note that <f) t (x,y) is orthogonal to 4> (x,y) by construction. Following this 
procedure we calculated the first four levels for the Henon-Heiles Hamiltonian. 
Our results are reported in Figs. 4-7. 
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We also calculated the variational ground state wave-function for this problem 
by minimizing the expectation value of the Hamiltonian, using an identical 
neural form. In Figs. 8-9, we plot the pointwise error, i.e. the (normalized) 
summand of eq. (11) for the collocation and the variational wavefunctions 
respectively. 

3. 6 Three Coupled Anharmonic Oscillators 

As a three-dimensional example we consider the potential for the three coupled 
sextic anharmonic oscillators [ |I~8|j : 

V(x, y, z) = V(x) + V(y) + V(z) + xy + xz + yz 

where 

V(x) = \x 2 + 2a; 4 + \x & 

The trial solution $ t {x,y,z) is parametrized as: 

<f> t (x, y, z) = e- x{x2+y2+z2) N{x, y, z, u, w {x \ w {y \w (z \ v), A > 

where N is a feedforward neural network with one hidden layer (with m = 25 
hidden units) and three input nodes (accepting the values of x, y and z): 

m 

N(x, y, z, u, w^ x \ w^ y \ w^ z \v) = ^2 VjCr(xWj X ^ + yWj + zWj + Uj) 

3=1 

We have considered a 28x28x28 grid in the [—4,4] x [—4,4] x [—4,4] domain 
both for computing the integrals and calculating the following error quantity 
that was minimized: 

/OO /"OO /"OO 
/ / (j) 2 t (x,y,z)dxdydz 
i.j.h- -oo J-oo J-oo 

where the energy is calculated by: 

= I-oo I-oo I-oo M x , V, z)H(f) t (x, y, z)dxdydz 
r» /-°°oo Ho y, z)dxdydz 

The ground state was computed and the corresponding eigenvalue was found 
equal to 2.9783, in agreement with the highly accurate result obtained by 
Kaluza |S . 
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4 Finite Element Approach 



The two-dimensional Schrondiger equation for the Henon-Heiles potential was 
also solved using the finite element approach in which the solution is expressed 
in terms of piecewise continuous biquadratic basis functions: 

^ = 5>i*i(£,n) (12) 
i=i 

where <3>j is the biquadratic basis function and ^ is the unknown at the i-th 
node of the element. The physical domain (x,y) is mapped on the computa- 
tional domain (£, n) through the isoparametric mapping: 

x = n) (13) 

i=i 

V = *EViMt,") (14) 

i=i 

where £ and n are the local coordinates in the computational domain (0 < 
£,n < 1) and Xj, the i-th node coordinates in the physical domain for the 
mapped element. 

The Galerkin Finite Element formulation calls for the weighted residuals Ri 
to vanish at each nodal point i — 1, . . . , N : 

Ri= I (Hip - ei>)§idet(3)didn = (15) 
Jn 

where J the Jacobian of the isoparametric mapping with 

a *fi\ dx d y dxd y 

det{3) = T^-^ (16) 

These requirements along with the imposed boundary conditions constitute a 
system of linear equations which can be written in a matrix form as: 

Kip = eMvp (17) 

where K is the stiffness and M is the mass matrix. The stiffness matrix in its 
local element form is: 

f [ f l r d$id&j <9<3>;<9$,, 1. 9 9 , ^ 1 , 9 1 ^ , , / ~w\ i j- j 

//¥rf + ^l^ + 2 ( * Wm^-/ -)*^}cfct(J)^ 

(18) 
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5x5 


7x7 


11x11 


16x16 


21x21 


29x29 


1.0075 


0.9997 


1.0015 


0.9994 


0.9989 


0.9986 


2.1988 


2.0852 


2.0037 


1.9930 


1.9911 


1.9901 


2.2001 


2.0862 


2.0037 


1.9930 


1.9911 


1.9901 


3.2495 


3.0159 


2.9767 


2.9648 


2.9593 


2.9571 


3.2878 


3.0515 


3.0065 


2.9943 


2.9885 


2.9857 


4.4347 


4.1139 


3.9868 


3.9433 


3.9323 


3.9262 



Table 1: Computed eigenvalues of the Henon-Heiles Hamiltonian using the 
FEM approach for various mesh sizes. 

The matrix M obtained above in its local element form is: 

QiQjdetffldtdn (19) 

Due to the Dirichlet boundary conditions zeros appear in the diagonal. Thus 
the mass matrix is singular and the total number of zeros in the diagonal of 
the global matrix is equal to the number of nodes on the boundaries and its 
degree of singularity depends on the size of the mesh. 

4-1 Extracting Eigenvalues and Eigenvectors 



/ 

Jn 



For the problem under discussion only the eigenvalues of the generalized eigen- 
value problem with the smallest real parts are needed. The eigenvalue problem 
is a symmetric one generalized eigenvalue problem but for generality purposes 
it is solved as a nonsymmetric one. Due to the size of the problem (from 1000 
- 4000 unknowns in our solution) direct methods are not suitable. 



We use Arnoldi's method as it has been implemented by Saad [[L4], [To], |16|j , 
which is based on an iterative deflated Arnoldi's algorithm. Saad proposes an 
iterative improvement of the eigenvectors as well as a Schur-Wiedland deflation 
to overcome cancellation errors in the orthonormalization of the eigenvectors 
at each step due to the finite arithmetic. 

If K is nonsingular, a simple way to handle the generalized eigenvalue problem 
is to consider the "reciprocal" problem: 

Mip = fjKip (20) 

where fi = 1/e. The infinite valued eigenvalues are transformed into zero 
eigenvalues. However, due to computer round-off errors, the infinite-valued 
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eigenvalues actually correspond to very large values in the calculations, which 
are turned into very small valued eigenvalues and not to exact zeros in the 
reciprocal problem. 

An alternative method would require the elimination of the rows with zero 
diagonal in the mass matrix, which are the rows corresponding essentialy to 
the boundary conditions. This scheme requires a number of manipulative op- 
erations on K and M which are prohibitive for large systems. The method is 
called the 'reduced algorithm' and requires the storage of the stiffness and the 
mass matrix. Other techniques have been proposed and mainly are transforma- 
tions of the generalized eigenvalue problem that map the infinite eigenvalues to 
one or more specified points in the complex plane [|17[]. The Shift-and-Invert 



transformation maps the infinite eigenvalues to zero. In the problem under 
discussion we have used the transformation: 

C = (M - <tK) _1 K (21) 

and the problem (20) is transformed to the problem: 

Ctp = fi'ifj (22) 

whose eigenvalues are related to those of equation (20) through the relation 
// = l/(fi — a), where a is a real number called shift. This transformation 
favors eigenvalues with real part close to the shift. The eigenvalues e of the 
original problem are then given by e = /// (1 + cr//). 

The generalized eigenvalue problem was solved on a rectangular domain. Figs. 
10 and 11 shows the evolution of the first and second eigenvalues as the num- 
ber of equidistributed elements of the mesh and consequently the number of 
unknowns increases. Convergence occurs for a grid of equal elements (29 x 29) 
which results in a system of 3481 unknowns. The convergence of the first six 
eigenvalues is also shown in Table 1. It is obvious that in order to get accurate 
eigenvalues, dense FEM meshes must be used and this limits the application 
of the method. 
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5 Conclusions 



We presented a novel method appropriate for solving eigenvalue problems of 
ordinary, partial and integrodifferential equations. We checked the accuracy 
of the method by comparing to a result that is analyticaly known, i.e. the 
ground state energy of the Morse Hamiltonian. We then applied the method 
to two realistic and interesting problems, namely to the Schrddinger and to 
the Dirac equations for a muonic atom. In these equations we take account of 
the finite protonic charge distribution as well as of the Vacuum Polarization 
effective potential. Preliminary calculations using a proton density delivered 
by Quasi- RPA have also been performed . Since both the Schrddinger and 
the Dirac equations can be solved analyticaly in the case of a point charge 
nucleus, (ignoring also the vacuum polarization correction), we conducted cal- 
culations (not reported in this article) and determined the energies for the 
4/ and 5g levels to within 1 ppm f20fl . The wide applicability of the method 
is shown by solving an integrodifferential problem, coming from the field of 
Nuclear Physics. The two dimensional benchmark, namely the Henon-Heiles 
Hamiltonian, that has been considered by many authors and solved by a host 
of methods, was considered as well. Here we obtained not only the ground 
state, but also some of the higher states, following a projection technique to 
supress the already calculated levels. Our results are in excellent agreement 
with the ones reported in the literature. We solved this problem also by a stan- 
dard finite element technique and we compared the computational resources 
and effort. It is clear that the present method is far more economical and effi- 
cient. Also, as we have previously shown [|TJ for the case of non-homogeneous 
equations, its interpolation capabilities are superb. Coming to an end, we 
solved a three-dimensional problem that imposes a heavier load. Again the 
results for the three-coupled anharmonic sextic oscillators are in agreement 
with the high precision ones obtained in [IS] by a semi-analytical method. 
The examples treated in this article are essentially single particle problems. 
(In example 3.4 the few-body nature is embeded in the non-local kernel). 
Many-body problems will impose a much heavier computational load, and 



hence the fast convergence property of the sigmoidal functions PHI , as well 
as the availability of specialized hardware become very important. Few-body 
problems may be handled by extending the method in a rather straightfor- 



15 



ward fashion. However for many-body problems it is not clear as of yet how 
to find a tractable neural form for the trial wavefunction. The method is new 
and of course there is room for further research and development. Issues that 
will occupy us in the future are optimal selection of the training set, networks 
with more than one hidden layers, radial basis function networks few-body 
systems and implementation on specialized neural hardware. 

We would like to acknowledge the anonymous referee for his useful suggestions 
that resulted in making the article more valuable. 
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Figure 1: Feedforward neural network with one hidden layer. 



Figure 2: Ground state of: (a),(b) the Dirac and (c) the Schrodinger equation 
for muonic atoms. 



Figure 3: Ground state of the non-local Schrodinger equation for the n + a 
system (e = -24.07644). 



Figure 4: Ground state of the Henon-Heiles problem (e = 0.99866). 



Figure 5: First excited state of the Henon-Heiles problem (e = 1.990107). 



Figure 6: Second excited state (degenerate) of the Henon-Heiles problem (e = 
1.990107). 



Figure 7: Third excited state of the Henon-Heiles problem (e = 2.957225). 



Figure 8: Pointwise normalized error for the collocation wavefunction. 



Figure 9: Pointwise normalized error for the variational wavefunction. 



Figure 10: Convergence of the first eigenvalue as a function of the mesh size 
(number of unknowns). 



Figure 11: Convergence of the second eigenvalue as a function of the mesh size 
(number of unknowns). 



18 




19 




20 




21 




22 




23 




24 




25 



0.0015 



0.001 



0.0005 







0.015 
0.01 
0.005 







1.04 - 




500 1000 1500 2000 2500 3000 



28 




29 



